In [1]:
import numpy as np
from spatio_flux import SpatioFluxVivarium, build_path, render_path

Make a Vivarium¶

view available types and processes

In [2]:
# TODO -- make a spatio-flux specific Vivarium with the the core preloaded
vi = SpatioFluxVivarium()
In [3]:
# view the available types
vi.get_types()
Out[3]:
['',
 'string',
 'length*mass/current^2*time^2',
 'mass/time^3',
 'list',
 'process',
 'length^2*mass/time^3',
 'reaction',
 'step',
 'bounds',
 'mass/current*time^2',
 'length^2*mass/temperature*time^2',
 'length^2*mass/time^2',
 'printing_unit',
 'length/time',
 'time',
 'length*mass/time^2',
 'current*time/substance',
 'length^4*mass/current*time^3',
 'emitter_mode',
 'number',
 'time^2/length',
 'kinetics',
 '/substance',
 'length^2',
 'current^2*time^4/length^2*mass',
 'path',
 'current*time',
 'current*length^2*time',
 'current*length^2',
 'temperature',
 'luminosity/length^2',
 'interval',
 'current^2*time^4/length^3*mass',
 'substrate_role',
 'length^0_5*mass^0_5/time',
 'length^1_5*mass^0_5/time',
 '/temperature*time',
 'tree',
 'wires',
 'length^2*mass/time',
 'quote',
 'current*length*time',
 'length/mass',
 'current*time/mass',
 'substance',
 'map',
 '/printing_unit',
 'tuple',
 'mass/length*time',
 'length^2*mass/current^2*time^2',
 '/time',
 'mass',
 'edge',
 'boundary_side',
 'current',
 'float',
 'length^3',
 'length^2*mass/current*time^3',
 'mass^0_5/length^1_5',
 'maybe',
 'enum',
 'positive_float',
 'length^3*mass/current^2*time^4',
 'length^2*mass/substance*temperature*time^2',
 'length*temperature',
 'length^1_5*mass^0_5/time^2',
 'luminosity',
 'length^2/time^2',
 'mass/temperature^4*time^3',
 'mass/length',
 'protocol',
 'time/length',
 'union',
 'schema',
 'current^2*time^3/length^2*mass',
 'boolean',
 'length*time/mass',
 'integer',
 'composite',
 'length',
 'length^2*mass/current^2*time^3',
 'length^4*mass/time^3',
 'array',
 'length^2/time',
 'substance/length^3',
 'printing_unit/length',
 '/length',
 'length^0_5*mass^0_5',
 'any',
 'length^3/mass*time^2',
 'current*time^2/length^2*mass',
 'mass^0_5/length^0_5*time',
 'length^3/time',
 'substance/time',
 'length^2*mass/current*time^2',
 'mass/length*time^2',
 'length*mass/current*time^3',
 'mass/length^3',
 'mass/time^2',
 'particle',
 'length/time^2']
In [4]:
# view the available processes
vi.get_processes()
Out[4]:
['composite',
 'ram-emitter',
 'DiffusionAdvection',
 'console-emitter',
 'DynamicFBA',
 'MinimalParticle',
 'Particles',
 'json-emitter']

dFBA¶

Dynamic Flux Balance Analysis (dFBA) extends traditional Flux Balance Analysis (FBA) to model the dynamic behavior of metabolic networks over time, allowing for the simulation of growth and substrate consumption in a changing environment.

In [5]:
# inspect the config schema for the 'increase' process
vi.process_config('DynamicFBA')
Out[5]:
{'model_file': 'string',
 'kinetic_params': 'map[tuple[float,float]]',
 'substrate_update_reactions': 'map[string]',
 'biomass_identifier': 'string',
 'bounds': 'map[bounds]'}
In [6]:
vi.core.default(vi.process_config('DynamicFBA')) # TODO use default to
Out[6]:
{'model_file': '',
 'kinetic_params': {},
 'substrate_update_reactions': {},
 'biomass_identifier': '',
 'bounds': {}}
In [7]:
# dfba_config  = vi.process_config('DynamicFBA', dataclass=True)  # TODO -- make this work
In [8]:
# TODO - enable get_dataclass to work with the new process
# dfba_config  = v.get_dataclass('DynamicFBA')
dfba_config = {
    "model_file": "textbook",
    "kinetic_params": {
        "glucose": (0.5, 1),
        "acetate": (0.5, 2)},
    "substrate_update_reactions": {
        "glucose": "EX_glc__D_e",
        "acetate": "EX_ac_e"},
    "biomass_identifier": "biomass",
    "bounds": {
        "EX_o2_e": {"lower": -2, "upper": None},
        "ATPM": {"lower": 1, "upper": 1}
    }
}
In [9]:
# make the vivarium
v = SpatioFluxVivarium()

# add a dynamic FBA process called 'dFBA'
v.add_process(name="dFBA",
              process_id="DynamicFBA",
              config=dfba_config,
              )
v.diagram(dpi='70')
Out[9]:
No description has been provided for this image
In [10]:
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]

for mol_id in mol_ids:
    v.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand()
    )

v.connect_process(
    name="dFBA",
    inputs={
            "substrates": {
                mol_id: build_path(path, mol_id)
                for mol_id in mol_ids}
        },
    outputs={
            "substrates": {
                mol_id: build_path(path, mol_id)
                for mol_id in mol_ids}
        }
)

# add an emitter to save results
v.add_emitter()

# TODO -- show_values does not work
v.diagram(dpi='70', show_values=True)
Out[10]:
No description has been provided for this image
In [11]:
v.set_value(path = ['fields', 'glucose'], value=10)
v.set_value(path = ['fields', 'biomass'], value=0.1)
field = v.get_value(['fields'])
print(field)

{'glucose': 10, 'acetate': 0.6199272885535164, 'biomass': 0.1}
In [12]:
v.save(filename='dFBA_t0')

Saved file: out/dFBA_t0.json
In [13]:
# run the simulation for 10 time units
v.run(interval=60)
In [14]:
v.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=[
        [
            '/fields/glucose',
            '/fields/acetate',
            '/fields/biomass'
        ]
    ]
)
Out[14]:
No description has been provided for this image

Spatial dFBA¶

In [15]:
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]
rows = 2
columns = 1

# make the vivarium
v2 = SpatioFluxVivarium()
for mol_id in mol_ids:
    v2.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand(rows, columns)
    )

# add a dynamic FBA process at every location
for i in range(rows):
    for j in range(columns):
        dfba_name = f"dFBA[{i},{j}]"

        # add a process for this location
        v2.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config,
        )
        v2.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}},
            outputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}}
        )

# add an emitter to save results
v2.add_emitter()

v2.diagram(dpi='70')
Out[15]:
No description has been provided for this image
In [16]:
# change some initial values
v2.merge_value(path = ['fields', 'glucose', 0, 0], value=10.0)
v2.merge_value(path = ['fields', 'biomass', 0, 0], value=0.1)
field = v2.get_value(['fields'])
print(field)

{'glucose': array([[10.        ],
       [ 0.13101474]]), 'acetate': array([[0.47702186],
       [0.37300945]]), 'biomass': array([[0.1       ],
       [0.56016474]])}
In [17]:
v2.run(60)
In [18]:
v2.get_timeseries(as_dataframe=True)
Out[18]:
/global_time /fields/glucose /fields/acetate /fields/biomass
0 0.0 [[10.0], [0.1310147415999353]] [[0.47702185787766227], [0.37300945027621435]] [[0.1], [0.5601647442497885]]
1 1.0 [[9.904761904761905], [0.014710265043333406]] [[0.4891299100737515], [0.0]] [[0.1079432568708625], [0.578804696512899]]
2 2.0 [[9.80200585245463], [0.014710265043333406]] [[0.5021057975498924], [0.0]] [[0.11651530697082768], [0.578804696512899]]
3 3.0 [[9.691145527078628], [0.014710265043333406]] [[0.5160006282993799], [0.0]] [[0.12576552148631978], [0.578804696512899]]
4 4.0 [[9.571550338512791], [0.014710265043333406]] [[0.530866005462576], [0.0]] [[0.13574706721093277], [0.578804696512899]]
... ... ... ... ...
56 56.0 [[0.0], [0.014710265043333406]] [[0.0], [0.0]] [[0.9748198839578434], [0.578804696512899]]
57 57.0 [[0.0], [0.014710265043333406]] [[0.0], [0.0]] [[0.9748198839578434], [0.578804696512899]]
58 58.0 [[0.0], [0.014710265043333406]] [[0.0], [0.0]] [[0.9748198839578434], [0.578804696512899]]
59 59.0 [[0.0], [0.014710265043333406]] [[0.0], [0.0]] [[0.9748198839578434], [0.578804696512899]]
60 60.0 [[0.0], [0.014710265043333406]] [[0.0], [0.0]] [[0.9748198839578434], [0.578804696512899]]

61 rows × 4 columns

In [19]:
# get a list of all the paths so they can be plotted together in a single graph
all_paths = []
for i in range(rows):
    for j in range(columns):
        # get the paths for this location
        location_path = []
        for mol_id in mol_ids:
            this_path = build_path(path, mol_id, i, j)
            rendered_path = render_path(this_path)
            location_path.append(rendered_path)
        all_paths.append(location_path)
# print(all_paths)
In [20]:
v2.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=all_paths
)
Out[20]:
No description has been provided for this image
In [21]:
v2.show_video()

Diffusion/Advection¶

This approach models the physical processes of diffusion and advection in two dimensions, providing a way to simulate how substances spread and are transported across a spatial domain, essential for understanding patterns of concentration over time and space.

In [22]:
vi.core.default(vi.process_config('DiffusionAdvection'))
Out[22]:
{'n_bins': (0, 0),
 'bounds': (0.0, 0.0),
 'default_diffusion_rate': 0.1,
 'default_diffusion_dt': 0.1,
 'diffusion_coeffs': {},
 'advection_coeffs': {}}
In [23]:
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
    'glucose',
    'acetate',
    'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
    'biomass': (0, -0.1)
}
path = ['fields']

v3 = SpatioFluxVivarium()

for mol_id in mol_ids:
    v3.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand(n_bins[0], n_bins[1])
    )

v3.add_process(name='diffusion_advection',
               process_id='DiffusionAdvection',
               config={
                   'n_bins': n_bins,
                   'bounds': bounds,
                   'default_diffusion_rate': diffusion_rate,
                   'default_diffusion_dt': diffusion_dt,
                   # 'diffusion_coeffs': diffusion_coeffs_all,
                   'advection_coeffs': advection_coeffs,
               },
               inputs={'fields': ['fields']},
               outputs={'fields': ['fields']}
               )

# add an emitter to save results
v3.add_emitter()

v3.diagram(dpi='70')
Out[23]:
No description has been provided for this image
In [24]:
v3.run(60)
In [25]:
v3.plot_snapshots(n_snapshots=5)
No description has been provided for this image
In [26]:
v3.show_video()

COMETS (Computation Of Microbial Ecosystems in Time and Space)¶

COMETS combines dynamic FBA with spatially resolved physical processes (like diffusion and advection) to simulate the growth, metabolism, and interaction of microbial communities within a structured two-dimensional environment, capturing both biological and physical complexities.

In [74]:
bounds = (20.0, 10.0)  # Bounds of the environment
n_bins = (12, 6)
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
    'biomass': (0, 0.1)
}
path = ['fields']

# Initialize the acetate concentration field with zeros
acetate_field = np.zeros((n_bins[0], n_bins[1]))

# Create a linear glucose concentration gradient from 1 (top) to 0 (bottom) vertically
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc

# Initialize the biomass concentration field, setting a small biomass concentration in a specific row
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1  # place some biomass

# make the vivarium
v4 = SpatioFluxVivarium()

# create the molecular fields
# for mol_id in mol_ids:
v4.add_object(
    name='glucose',
    path=path,
    value=glc_field
)
v4.add_object(
    name='biomass',
    path=path,
    value=biomass_field
)
v4.add_object(
    name='acetate',
    path=path,
    value=acetate_field
)

# add a diffusion/advection process
v4.add_process(name='diffusion_advection',
               process_id='DiffusionAdvection',
               config={
                   'n_bins': n_bins,
                   'bounds': bounds,
                   'default_diffusion_rate': diffusion_rate,
                   'default_diffusion_dt': diffusion_dt,
                   'advection_coeffs': advection_coeffs},
               inputs={'fields': path},
               outputs={'fields': path}
               )

# add a dynamic FBA process at every location
for i in range(n_bins[0]):
    for j in range(n_bins[1]):
        dfba_name = f"dFBA[{i},{j}]"

        # add a process for this location
        v4.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config,
        )
        v4.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}},
            outputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}}
        )

# add an emitter to save results
v4.add_emitter()
In [75]:
v4.diagram(dpi='70',
           remove_nodes=[f"/dFBA[{i},{j}]" for i in range(n_bins[0]-1) for j in range(n_bins[1])]
           )
Out[75]:
No description has been provided for this image
In [76]:
v4.run(60)
In [77]:
v4.plot_snapshots(n_snapshots=5)
No description has been provided for this image
In [78]:
v4.show_video()
In [79]:
v4.plot_timeseries(
    subplot_size=(8, 3),
    query=[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0',
    ],
    combined_vars=[[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0',
    ]]
)
Out[79]:
No description has been provided for this image

Particles¶

In [80]:
# TODO -- remove this from final
import numpy as np
from spatio_flux import SpatioFluxVivarium
In [81]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (20, 40)  # Number of bins in the x and y directions

v5 = SpatioFluxVivarium()

v5.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.4,
        'boundary_to_add': ['top']
    },
)
v5.connect_process(
    name='particle_movement',
    inputs={
        'fields': ['fields'],
        'particles': ['particles']
    },
    outputs={
        'fields': ['fields'],
        'particles': ['particles']
    }
)

v5.initialize_process(
    path='particle_movement',
    config={'n_particles': 2}
)

v5.add_emitter()
v5.diagram(dpi='70')
Out[81]:
No description has been provided for this image
In [82]:
v5.run(100)
v5_results = v5.get_results()
In [83]:
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif

# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
    v5_results,
    skip_frames=3,
    bounds=bounds
)

Saving GIF to species_distribution_with_particles.gif
In [84]:
v5.diagram(dpi='70')
Out[84]:
No description has been provided for this image
In [85]:
v5.save(filename='v5_post_run', outdir='out')

Saved file: out/v5_post_run.json

Minimal particle process¶

In [86]:
from spatio_flux.processes.particles import get_minimal_particle_composition

particle_schema = get_minimal_particle_composition(v5.core)

document = v5.make_document()
# document['state']['particles'] = {}
# # document['composition']['particles']
# document['composition'] = particle_schema
document['state']['global_time'] = 0.0

v6 = SpatioFluxVivarium(document=document)

v6.diagram(dpi='70')
Out[86]:
No description has been provided for this image
In [87]:
v6.run(10)
v6_results = v6.get_results()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[87], line 1
----> 1 v6.run(10)
      2 v6_results = v6.get_results()

File ~/code/vivarium-interface/vivarium/vivarium.py:564, in Vivarium.run(self, interval)
    561 if not self.emitter_paths:
    562     self.add_emitter()
--> 564 self.composite.run(interval)

File ~/code/process-bigraph/process_bigraph/composite.py:783, in Composite.run(self, interval, force_complete)
    781 for path in self.process_paths:
    782     process = get_path(self.state, path)
--> 783     full_step = self.run_process(
    784         path,
    785         process,
    786         end_time,
    787         full_step,
    788         force_complete)
    790 # apply updates based on process times in self.front
    791 if full_step == math.inf:
    792     # no processes ran, jump to next process

File ~/code/process-bigraph/process_bigraph/composite.py:680, in Composite.run_process(self, path, process, end_time, full_step, force_complete)
    677     full_step = interval
    679 if future <= end_time:
--> 680     update = self.process_update(
    681         path,
    682         process,
    683         state,
    684         process_interval
    685     )
    687     # update front, to be applied at its projected time
    688     self.front[path]['time'] = future

File ~/code/process-bigraph/process_bigraph/composite.py:622, in Composite.process_update(self, path, process, states, interval, ports_key)
    600 """Start generating a process's update.
    601 
    602 This function is similar to :py:meth:`_invoke_process` except in
   (...)
    617     ``store``.
    618 """
    620 states = strip_schema_keys(states)
--> 622 update = process['instance'].invoke(states, interval)
    624 def defer_project(update, args):
    625     schema, state, path = args

File ~/code/process-bigraph/process_bigraph/composite.py:303, in Process.invoke(self, state, interval)
    302 def invoke(self, state, interval):
--> 303     update = self.update(state, interval)
    304     sync = SyncUpdate(update)
    305     return sync

File ~/code/spatio-flux/spatio_flux/processes/particles.py:183, in Particles.update(self, state, interval)
    180 # Apply diffusion and advection
    181 dx, dy = np.random.normal(0, self.config['diffusion_rate'], 2) + self.config['advection_rate']
--> 183 new_x_position = particle['position'][0] + dx
    184 new_y_position = particle['position'][1] + dy
    186 # Check and remove particles if they hit specified boundaries

TypeError: can only concatenate str (not "numpy.float64") to str
In [ ]:
from spatio_flux.processes.particles import get_minimal_particle_composition

bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (4, 8)  # Number of bins in the x and y directions

# same as before
v6 = SpatioFluxVivarium()
v6.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.4,
        'boundary_to_add': ['top']})
v6.connect_process(
    name='particle_movement',
    inputs={
        'fields': ['fields'],
        'particles': ['particles']},
    outputs={
        'fields': ['fields'],
        'particles': ['particles']})

# Here we go beyond the previous particle simulation
# add a minimal particle process into the schema
particle_schema = get_minimal_particle_composition(v6.core)
v6.set_schema('particles', particle_schema['particles'])

v6.initialize_process(
    path='particle_movement',
    config={'n_particles': 1})

v6.diagram(dpi='70')
In [ ]:
# v6
In [ ]:
# particle_schema = get_minimal_particle_composition(v6.core)
# particle_schema
In [ ]:
v6.run(10)
v6_results = v6.get_results()
In [ ]:
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif

# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
    v6_results,
    skip_frames=2,
    bounds=bounds
)